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SUMMARY 


An accurate shock capturing scheme Is developed for high speed flow. The 
flux limited dissipation scheme can be converted to a symmetric TVD scheme In 
the case of a scalar conservation law. Application to an oblique shock wave 
problem shows the Improved resolution and accuracy when compared to the 
results obtained by the adaptive dissipation method, while exhibiting com- 
parable computational efficiency. 


INTRODUCTION 


It has been a difficult task for computational aerodynamiclsts to capture 
the nonoscl 1 latory shock waves. Unbounded growth of spurious oscillations 
often resulted In numerical Instability. Central difference schemes produce 
good solutions In smooth regions of the flow, but are prone to oscillations In 
the neighborhood of shock waves. It Is well-known that these wiggles can be 
suppressed by the Introduction of a carefully controlled blend of first and 
third order numerical dissipation terms. Recently the authors showed that the 
strong oblique shock waves In high Mach number flows can be captured by the 
method of adaptive dissipation (ref. 1). 

It has long been recognized that upwind differencing can eliminate unde- 
sirable oscillations near shock waves. Stemming from the mathematical theory 
of scalar conservation laws, Harten proposed the concept of total variation 
diminishing (TVD) schemes (ref. 2). TVD schemes preserve the monotonicity of 
an Initially monotone profile, because the total variation would Increase If 
the profile ceased to be monotone. Second order schemes can be constructed by 
using multipoint extrapolation formulas to estimate the numerical flux, or by 
adding higher order dissipative terms. In either case flux limiters are then 
needed to control the signs of the coefficients of a semi-discrete approxima- 
tion to the hyperbolic system (refs. 2 to 5). However, TVD schemes In two or 
more space dimensions are at most first order accurate (ref. 6). Nevertheless, 
two-dimensional results using one-dimensional second order TVD schemes and 
dimensional splitting show sharp resolution of discontinuities without oscil- 
lations (ref. 7). Jameson constructed an efficient scheme which Is converted 
to a TVD scheme In the case of a scalar conservation law using flux limited 
dissipation (ref. 5). In this paper, we Introduce a modified form of the 
scheme and apply It to high Mach number flow calculations. Computed solutions 
are compared with the results obtained by the adaptive dissipation. The LU 
Implicit scheme (ref. 1) Is used as a baseline algorithm. 


GOVERNING EQUATIONS 


The Euler equations are obtained from the Navler-Stokes equations by 
neglecting viscous terms. Let p, u, v, E, H, and p be the density, Car- 
tesian velocity components, total energy, total enthalpy, and pressure, and let 
x and y be Cartesian coordinates. Then for a two-dimensional flow these 
equations can be written as 


aw 3 F 3G 
at + ax + ay 


0 


( 1 ) 


where W is the vector of dependent variables, and F and G are convective 
flux vectors 


W = (p,pu,pv,pE) T 

2 T 

F = ( pu , pu + p,pvu,u(pE + p)) (2) 

2 T 

G = ( pv , puv ,pv * p,v(pE * p)) 

The pressure Is obtained from the equation of state 

P = p(y - 1){e - \ (u 2 * v 2 )| (3) 

These equations are to be solved for a steady state aw/at = 0 
where t denotes time. 


SEMI-DISCRETE FINITE VOLUME METHOD 


A convenient way to assure a steady state solution Independent of the time 
step Is to separate the space and time discretization procedures. In semi- 
discrete finite volume method one begins by applying a semi-discretization in 
which only the spatial derivatives are approximated. The use of a finite vol- 
ume method for space discretization allows one to handle arbitrary geometries 
and helps one to avoid problems with metric singularities that are usually 
associated with finite difference methods. The scheme reduces to a central 
difference scheme on a Cartesian grid, and Is second order accurate In space 
provided that the mesh Is smooth enough. It also has the property that uniform 
flow Is an exact solution of the difference equations. 


LU IMPLICIT SCHEME 


Let the Jacobian matrices be 



B 


aG 

aw 


2 


and let the correction be 


6W - W 


n+1 


W 


where n denotes the time level. 


The linearized Implicit scheme for a system of nonlinear hyperbolic equa- 
tions such as the Euler equations can be formulated as . 


{i * B At (D x A * D y B)j 


6W + At R = 0 


(4) 


where R Is the residual 


R = D x F(W n ) + D y G(w") 

Here D x and Dy are central difference operators that approximate 3/3x 
and 3/dy. 

If B = 1/2 the scheme remains second order accurate In time, for other 
values of 8, the time accuracy drops to first order. The unfactored Implicit 
scheme (eq. (4)) produces a large block banded matrix which Is very costly to 
Invert and requires huge storage. An unconditionally stable Implicit scheme 
which has error terms at most of order (At) 2 In any number of space dimen- 
sions can be derived by the LU factorization. 

I * 8 At (D'A + + d'bM Jl + 8 At (oV + dV) 
x y | / * y 

6W + At R = 0 (5) 


where D x and Dy are backward difference operators and D x 

and DJ are forward difference operators. The reason for splitting Is 

to ensure the diagonal dominance of lower and upper factors as well as to make 

use of the built-in Implicit dissipation. 

Here, A + , A - , B + , and B are constructed so that the eigenvalues of " + " 
matrices are nonnegative and those of matrices are nonpositive. 


where 


A + = ^ (A f r A I), A - = ^ (A - r A I) 

B + = 2 (B * r B I), B~ = ~ (B - r B I) 


( 6 ) 


r A > max( 1^1). 


r B ^ max( | x g | ) 


(7) 


Here, and Xg represent eigenvalues of Jacobian matrices. Equation (5) 
can be Inverted In two steps. The LU Implicit scheme needs the Inversion of 
sparse triangular matrices which can be done efficiently without using large 
storage. This scheme has only two factors In three dimensions. 
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FLUX LIMITED DISSIPATION SCHEME 


Consider the scalar conservation law 


it + k f ( u > =° 


( 8 ) 


It Is well known that the total variation 


TV . / I I d* 


cannot Increase. Suppose now that a multi-point semi-discrete approximation 
to equation (8) is expressed in the form 


du Q-l 

T7 = C ( 1 ) ( U. - U. -. ) 

dt qT-q q M M-V 


(9) 


The discrete total variation is 


TV . E lu 1 - u, 

1 — — 00 

It can be shown (ref. 8) that this will not Increase If and only if 


and 


c^d-l) > c _ 2 0-2). . . > c Q ( 1 -Q) > 0 
-C 0 (1) > - c ^ ( 1 +1 ) . . . > -c (HQ-1) > 0 


( 10 ) 


Now consider the scheme 

du 


dt + Ax (h Hl/2 " h 1-l/2 ) “ 0 


( 11 ) 


where h^ +1 ^ 2 H an approximation to the flux across the boundary between the 
(Hl) s * and l*^ 1 cells. Denoting f(u^) by f^, define the numerical flux as 


h Hl/2 = 2 (f Hl * f 1 ) + d Hl/2 


( 12 ) 


where d^ 2 a d Hs1pative flux. 


d 1+l/2 a 1+l/2 1 B(e i+3/2’ e i+l/2 ) ' 2 e i+l/2 + B(e i+l/2’ ®1-1/2N 


(13) 


where e. . .„ = p. , - p. for example. 
Hl/2 1+1 1 


B (p,q) is known as Roe's min mod function (ref. 4) and can be defined as 

B(p,q) =js(p) + s(q) j min ( |P| , |q|) (14) 


V 



where 


Here 


s(p) 


| If P i 0 
If P < 0 


‘ 1 + 1/2 


- min (1/2, k 0 + k 1 l V-|/ 2 )( R i +1 R ^) 


(15) 


where k and k, are the constants and 
o 1 


1+1/2 

- max(v ul/2 , 

w 1+l 


P, 

2P. 

P, , 


1+1 - 

1 + 

1-1 

V 1 = 

P 1+l * 2P 1 + 

P 1-l 


The spectral radius R can be estimated as the value of 

R = | ay u - axv | + cyax 2 + ay 2 

on the edge separating cells (1+1, j) and (1,j) In two-dimensions, c Is the 
speed of sound. 


RESULTS 


A two-dimensional schematic high speed Inlet Is selected as a model prob- 
lem to test the performance of the new shock capturing scheme. A 54 by 32 
H-mesh Is shown In figure 1. Figure 2 shows the plot of Mach number contours 
for the Mach 5 flow using the adaptive dissipation method. Improved result by 
the TVD scheme Is shown In figure 3. The plots of Mach number along the cen- 
terline are compared In figure 4 for the adaptive dissipation method and In 
figure 5 for the TVD scheme. While both methods capture nonoscl 1 latory shock 
waves, the results of the TVD scheme show sharper resolution and Improved 
accuracy. Considering that the computational mesh lines are not aligned with 
the oblique shock wave, the result seems to be satisfactory. The convergence 
rate for the TVD scheme (fig. 7) Is found to be comparable to that for the 
adaptive dissipation method (fig. 6). Moreover, the computational work per 
cycle of the present symmetric TVD scheme Is comparable to that of the adaptive 
dissipation scheme. In general, symmetric TVD schemes offer a significant 
reduction of computational complexity In comparlsm with upwind TVD schemes, 
while exhibiting comparable shock capturing capabilities. 
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Figure 3. - Mach number contours for Mach 5 flow (present method) 
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Figure - Mach number along the centerline (adaptive dissipation 

SCHEME). 
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